Lattice Boltzmann-Langevin Equations 11 

J.W. Duftytand M.H. Ernst* 
February 6, 2008 



Abstract 

Intrinsic fluctuations around the solution of the lattice Boltzmann 
equation are described or modeled by addition of a white Gaussian noise 
source. For stationary states a fluctuation-dissipation theorem relates the 
variance of the fluctuations to the linearized Boltzmann collision operator 
and the pair correlation function. 



1 Introduction 

The discussions and contributions at this meeting have shown a growing interest 
in fluctuations both in the microscopic lattice gas approach, as well as in the 
phenomcnological lattice Boltzmann approach. The emphasis was on unstable 
systems (phase separation, spinodal decomposition, pattern formation) where 
fluctuations drive the system away from a spatially uniform equilibrium state. 
So far the noise considered in connection with the lattice Boltzmann approach 
was introduced in a rather ad hoc fashion. It may be additive or multiplicative, 
external or intrinsic. There do not seem to be any systematic studies that make 
the connection with the multitude of results for fluctuations and transport in 
the continuous case (see Bixon and Zwanzig [1969], Fox and Uhlenbeck [1970], 
Ernst and Cohen [1981] and Marchetti and Dufty [1983]). The goal of this paper 
is to provide this link and to show that the well-etablished results on intrinsic 
fluctuations, derived from the Boltzmann-Langevin equation in the continuous 
case, can be extended to the lattice Boltzmann equation. 

The Boltzmann equation is a deterministic equation for the average occu- 
pation f(r, c, t) of a one particle state specified by position r and velocity c. 
When derived from an underlying microdynamics, it describes the 'slow' time 
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evolution of / (r, c, t) on coarse grained spatial and temporal scales of the or- 
der of the mean free path £ and time t between collisions. The deviations 
of the mesoscopic occupation, n(r,c,t), from this average Boltzmann value re- 
flect fluctuations due to other 'fast' degrees of freedom that are averaged out 
in obtaining the Boltzmann equation. A measure of these fluctuations is given 
by their correlation functions, which also can be calculated from the micrody- 
namics. A unified picture of both the transport and fluctuations is given by a 
Boltzmann-Langevin (B-L) equation. The existence of other degrees of freedom 
is recognized in the B-L equation by adding a stochastic source that generates 
fluctuations in the solutions to this equation relative to the deterministic so- 
lution. This is a particularly convenient representation for simulation of noise 
effects on the Boltzmann equation. The characteristics of this stochastic source 
should be determined from the microdynamics for consistency with the direct 
evaluation of correlation functions using methods of statistical mechanics and 
kinetic theory. Studies along these lines are in progress [Dufty and Ernst 1994]. 

In other contexts, a Boltzmann equation is simply postulated as a mathemat- 
ical model to simulate and study macroscopic fluid properties. Its extension to a 
corresponding B-L equation to describe the effects of noise now has no underly- 
ing microdynamics to characterize that noise. To assist in this characterization, 
we establish in the next section the general relationship of the noise variance to 
the pair correlation function at equal time and different times. The form of this 
relationship is independent of any underlying microdynamics. This connection 
of the noise to correlation functions provides a somewhat more physical basis 
for suggesting the noise characteristics. Furthermore, known results about the 
correlation functions for systems with a microdynamics can be 'borrowed' for 
application to the mathematical modeling as well. 

Our primary interest here is to make the connections just mentioned and to 
discuss their content at a phenomenological level. The B-L equation is defined 
in the next section and it is assumed that the noise is white (uncorrelated 
in time) and essentially characterized by its variance (this is strictly the case 
only for Gaussian noise). Equations for the correlation of fluctuations at equal 
and different times are derived from the B-L equation, under the condition 
that the fluctuations are small. It is observed that these equations are the 
precise analog of those derived for a real low density gas (see Ernst and Cohen 
[1981], Marchetti and Dufty [1983]). For the special case of fluctuations about a 
stationary state a fluctuation-dissipation relation is obtained, defining the noise 
variance in terms of the stationary state correlations and the linear Boltzmann 
operator. For non-stationary states additional information is required. The 
analysis applies as well to unstable states up to times for which the fluctuations 
remain small. In section 2 we give a concrete example of a BGK-Langevin 
description of fluctuations about a steady state. Finally, in section 3 we show 
that hydrodynamic equations with noise can be obtained consistently from the 
B-L description. The expected Landau-Lifshitz form is confirmed. 
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2 B-L Equation and Correlation Functions 



In the absence of noise, the lattice Boltzmann equation describes deterministic 
dynamics for the occupation number, n^r, t), of a single particle state at node 
r in velocity channel Cj where i labels the allowed velocity states. The rii(r,t) 
are continuous real variables on the interval [0, a] where a is the maximum 
occupancy of a state (in the case of an underlying microdynamics with Fermi 
exclusion rules, a value a > 1 implies that these occupation numbers are coarse 
grained averages over several nodes). The corresponding B-L equation with 
noise is, 

n(x, t + 1) = V(x | n{t)) + S x F(x, t), (1) 

where a simplified notation n(x, t) = rii(r, t) has been introduced. The generator 
for the deterministic dynamics, V(x \ n(t)), is given by 

V(x | n(t)) = S x (n(x, t) + I(x \ n(t))), (2) 

where I(x | n(tj) — Ii(r \ n{t)) is the non-linear Boltzmann collision operator, 
and S x is the free streaming operator defined by S x f(f) = f(r — Cj) for any 
function of r. Finally, the stochastic source, F(x,t), is assumed to represent 
Gaussian white noise with zero mean value, statistically independent of the 
occupation numbers, (n(x,t)F(y,T)) — for t < t, and covariance specified by 
a noise intensity matrix, B(x,y), 

(F(x, t)F(y, t 1 )) = S(t, t')B(x, y | n{t)). (3) 

The average is taken over the noise source. Furthermore the noise does not 
generate any mass, momentum or energy, so that F is orthogonal to the sum- 
mational invariants. We note that the noise matrix B(x,y \ n) in general can 
depend on the state at time t; this is an example of multiplicative noise. Once 
the noise matrix has been specified, the above provides a complete description 
of fluctuations and transport. 

To simplify the notation, let n denote the vector whose components are 
n(x,t). Similarly, V(n), F, and B(n) denote the vectors and matrix with 
elements V(x\n), F(x), and B(x,y \ n), respectively. The B-L equation is then, 

n(t + 1) = V(n(t)) + SF(t), (F(t)F(f )) = 8(t, f )B(n(t)). (4) 

The solution to the deterministic (noise- free) Boltzmann equation, f(t), is ob- 
tained from 

f(t+l)=V(f(t)). (5) 

The fluctuations around this deterministic solution are measured by <5n(t) = 
n(t) — f(t). The corresponding pair correlation functions are defined by, 

g(t) = (Sn(t)5n(t)), c(t,t') = (Sn(t)Sn(t')). (6) 
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Closed equations for these correlation functions follow in two steps: first the 
B-L equation is applied to g(t + 1) and c(t + 1, t'); next the resulting averages 
are expanded to lowest order in 6n(t), for small fluctuations, 



g(t + l) = L(t)g(£)L T (t) + SB(i)S T 



(7) 



c(t + l,t') = L(t)c(t,t'). 
Here B(t) = B(f (t)) and L(t) is the matrix, 



(8) 



L(i) = -V(n) | f(t) . 



(9) 



Both the equal time and two time correlation functions obey linear equations 
governed by the matrix L(t) whose form is determined from the generator, V(f), 
for the given Boltzmann equation. They depend on the specific nonequilibrium 
solution to the Boltzmann equation, f(t), being considered. For the special case 
where f (i) = f is the equilibrium state, this result forms the starting point 
in the calculation of the Green-Kubo formulas for lattice gas cellular automata 
(LGCA) in the mean field approximation [Brito et al 1991], as well as for the 
dynamic structure factor in such systems [Grosfils et al 1993]. 

We remark at this point that ([s]) and ([?]) through are the direct analog of 
corresponding results from the kinetic theory of real gases at low density. The 
difference here is that the form of B(t) is given explicitly from the theory there, 
whereas here it is as yet unspecified. 

Before discussing the content of these equations it may be useful to see how 
they can be applied to determine the noise matrix in a special case. Consider 
fluctuations around a stationary (or quasi-stationary) solution to the Boltzmann 
equation, f . For consistency, the noise matrix must be chosen such that the 
equal time correlation functions are also stationary. Consequently, (Q) gives, 



SB S T = g - L og L^ or B = S T gS - (1 + H)g(l + fi T ), (10) 



where B Q and L D are given by B and L evaluated at the stationary state f Q , 
and L Q = S(l + fl) has been expressed in terms of the linearized collision 
operator fl. This is the classical fluctuation-dissipation theorem extended to 
lattice gas automata. It expresses the noise matrix in terms of the stationary 
state correlations and the Boltzmann collision operator linearized around the 
stationary state. A more explicit form is given in the example at the end of this 



How would ( [To|) be used to implement a simulation of noise in a B-L equa- 
tion? First the model for the Boltzmann equation is specified by using collision 
rules for a LGCA or by a BGK model; next the equal time correlation matrix is 
determined for the model; finally, solutions to the B-L equation are simulated 
for initial states near f Q , selecting the noise F from a Gaussian whose half width 
B Q is determined from the fluctuation-dissipation relation. In this simulation of 



section. 
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the B-L equation in the stationary state the equal time correlation function g is 
used as input to determine the noise statistics B , and the two time correlation 
function c(t, t') = Sn(t)Sn(t') can be measured by simulations, even for models 
without any underlying microdynamics. The simulated values should agree with 
the theoretical value calculated from (||). 

To illustrate the procedure sketched below (|l(]) we consider two examples: a 
standard LGCA that obeys the conditions of semi-detailed balance (SDB), and a 
BGK-model. In the LGCA all equilibrium fluctuations are totally uncorrelated, 
i.e., 

g(x,y) = S(x,y)gi (x = f i;y = f'j) (11) 

with a known value, gi = fo%{l — foi) for the fluctuations in a single particle state 
because of the Fermi exclusion rule. Furthermore, the collision rules are strictly 
local such that the collision matrix in ([ll]) has the form Sl(x,y) = <5(r, r')tlij. 
The noise strength (10) reduces to 



B(x,y) = 8(r,r')B ij = -<5(r, r') [Sl^gj + + n ie g e fl je } 



(12) 



Let us compare these results with those of Bixon and Zwanzig [1969], and Fox 
and Uhlenbeck [1970] for the strength of the fluctuating term in the continuous 
case. The first two terms are direct analogs of the corresponding terms in the 
continuous case. The term quadratic in O is a typical lattice effect, similar to the 
so called 'propagation part' of the transport coefficients, as occurring in lattice 
gas automata (see e.g. Ernst and Dufty [1991]) and in the lattice Boltzmann 
approach. This will be shown in section 3. 

As a second example we consider a typical athermal BGK-model equation, 
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fi(f + Ci,t + 1) - fi(r, t) = — [Mr, t) - //] , 



(13) 



where the local equilibrium distribution ff depends on the local density p(r, t) = 
Si/i(f, t) and the local flow velocity p(r,t)u(r,t) — Y. i c i f i (f,t). We linearize 
this equation around basic equilibrium f i with vanishing flow velocity, where 
foi equals f a for a rest particle (i = 0) and equals / for a moving particle. The 
collision operator is diagonal and can be identified as, 



dp p du 



(14) 



where the subscript (o) indicates that the derivatives are taken at u = 0. 

As a specific example we take the BGK-model defined in Eq.(15) of Ernst 
and Bussemaker (see these proceedings) and obtain, 



d 



(cl + cl) 



(15) 
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where i refers to a moving particle. Furthermore, d is the spatial dimension of 
the lattice, | &i |= c G the lattice distance, and b the coordination number. We 
also used the speed of sound from the above reference, i.e., 

c 2 = dp = bc ° 2 ^ = ° 2 ° (1 d f°) (16) 
5 dp d dp d dp 

So far we have specified the matrix elements of the collision operator to be used 
in the expression (|l2|) for the noise strength. Next we address the choice of 
equilibrium fluctuation <?.; = ((Srii) 2 ). In lattice gas automata the exclusion rule 
for occupation numbers makes the equilibrium fluctuations gi similar to those 
for an ideal Fermi gas. In the mathematical modeling of the lattice Boltzmann 
approach the positive function gi can be chosen freely by lack of an underlying 
microscopic model. However, the choice gi = p(df°/dp), which resembles the 
result for a Maxwell-Boltzmann gas, is closest in spirit to the BGK-model. This 
is so because the local equilibrium distribution function //, used in typical lattice 
gas applications [Chen et al 1992], depends on the local average density p and 
the fluid flow velocity u, and resembles a local Maxwell-Boltzmann distribution 
expanded to terms of 0{u 2 ) for small velocities. In the typical model under 
consideration we have for the rest particles g = p(l — dc 2 s /c 2 D ) and for moving 
particles gi — (pd/b){c 2 /c 2 ). This completes the determination of the linear 
BGK-Langevin equation. 



3 Noise in Fluid Dynamics 

In this section we calculate the correlation function in equilibrium of the fluc- 
tuating part of the single particle distribution function in the limit of large 
spatial and temporal scales, for systems where the noise strength is determined 
by B (x,y) in ( 12]) . From that result one can compute the fluctuations in the 
pressure tensor or in the heat current, if energy is also conserved in the lattice 
Boltzmann equation. The results obtained will be compared with the Landau- 
Lifshitz formulas for the correlation strength in fluctuating hydrodynamics. 
For the analysis of this section it is convenient to work with Fourier and Laplace 
transforms, defined through, 

oo oo 

hi(k,z) =Y^Y. e ~ lt? ~ Zt6n ^t) =5>" z *n 4 (M). (17) 

t=0 r t=0 

For fluctuations around equilibrium the B-L equation becomes, 

e^+^fc, z) = (% + %)ni(fc, z) + Fi{k, z). (18) 
In vector and matrix notation the formal solution of ( |l8| ) reads 

n(k,z) = n (k,z) + f(k,z), (19) 
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where the sure (deterministic) part is ho(k,z) = A(k, z) e lk ' c+z n{k, 0) and the 
fluctuating part is f(k,z) = A(k, z) F(k, z). Here Aij(k,z) and [exp(ifc ■ c)]ij 
= Sij exp(ik • ci) are matrices, and h(k,z), etc are vectors with components 
fii(k, z), etc. The matrix A is defined as 

A(k,z) = [ e its+z - 1 -ftp 1 . (20) 

The sure part of (|l^) approaches the Chapman-Enskog hydrodynamic solution 
of the Boltzmann equation for large spatial and temporal scales. The fluctuating 
part of the distribution function has correlations given by, 

V~ 1 {fi(k, z)f*(k, z')) = V- 1 A u (k,z){F e (k,z)F^(k,z'))A* m (k,z'), (21) 

where the asterisk represent complex conjugation and where V is the number of 
nodes in the lattice. The correlation function of the F's is obtained by taking 
Fourier and Laplace transforms of (g), and using translational invariance with 
the result, 

V- 1 (fi(k,z)j;(k,z')) = Au(k,z)B lm A* jm (k,z')/[l-e- z - z ']. (22) 

As we are interested in hydrodynamic space and time scales, we take the limit 
k — > 0,z — > 0. Then A(k,z) — > — I/O, where the inverse of ft is only defined 
in the subspace orthogonal to the zero-eigenfunctions, the so-called collisional 
invariants. Next we invert the transforms in ( p2] ) to find the fluctuations fi (r, t) 
in the distribution function as, 

(Mr,t)f*(f',t')} ~ i(f,f')i(M')^^^ = 
-6(r,r>)5(t,t>) + \).. 9j + (i + \) .. 9i ] . (23) 



The second equality is a consequence of (p_2|). This is the final result for the fluc- 
tuations in the distribution functions on hydrodynamic space and time scales. 
The equation also shows that the term oc ft 2 in ( |l0| ) and ( |l2| ) directly trans- 
forms into the two terms in ( |23| ) containing |, which constitute the so called 
'propagation' part of the transport coefficients. 



The result ([23|) enables us to calculate the strength of the fluctuations in 
the dissipative currents in the fluctuating hydrodynamic equations, such as the 
stress tensor or the heat current in lattice gases with energy conservation. We 
illustrate this by calculating the fluctuations in the shear stress tl xy (r , t) — 

c X iC y ifi(r,t) with the help of (|23|). This yields, 

(IL xy (r,t)fl xy {r\t'))=2(c x \c x )v5(r,r')6(t,t'), (24) 
where v is the kinematic viscosity, given by 

v = -{c x c y \ f^ + ^J I c x c y ) / (c x \c x ) . (25) 
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The brackets are defined as, 

v= (a | M \b) =a i M ij g j b j , (26) 

which includes the definition of the inner product for M = 1. This expression 
is the standard formula for the kinematic viscosity as obtained from a lattice 
Boltzmann equation [Brito et al 1991]. In lattice gas automata the occupation 
numbers obey an exclusion principle, sothat g t = {{5ni) 2 ) — f°(l — f°). In 
athermal lattice gases with or without rest particles the single particle state 
distribution function f° — f is independent of the speed \ci\ and equal to the 
average density p per node, divided by the number of allowed velocity channels. 

4 Concluding Remarks 

We offer several remarks to summarize and put the above results in context. 

1) The noise considered here has been assumed to be white and Gaussian, 
and the resulting Boltzmann-Langevin equations define a discrete time Markov 
process. It is entirely characterized once the deterministic dynamics (i.e., the 
form of the Boltzmann collision operator V) and the noise intensity B is spec- 
ified. Non-Gaussian noise could be considered as well, and most of the results 
obtained here for pair correlations still hold. However, higher order correlations 
would be different. 

2) For stationary and quasi-stationary states the noise intensity has been 
determined by the consistency between the solution to the linearized Boltzmann- 
Langevin equation for equal time pair correlations and a specified stationary 
value for this correlation. The important new result here is the fluctuation- 
dissipation theorem (jl^). It relates the fluctuations to the matrix elements 
of the linear collision operator and the equilibrium value of the correlations 
g{x 1 y), which are assumed to be given. The relation is somewhat similar to that 
for continuous gases, but there are some typical lattice effects, caused by the 
discreteness of time and closely related to the so called 'propagation' viscosities 
of lattice gases. 

3) The linear Boltzmann-Langevin equations imply associated hydrodynamic 
equations with noise. The latter occurs as an additional component to the stress 
tensor. The fluctuation formula for the stress tensor (pi[), as derived from the 
Langevin noise added to the lattice Boltzmann equation, is in complete agree- 
ment with the results of Landau and Lifschitz for the noise in fluid dynamics. 
Furthermore, the explicit form of (^4|) identifies the noise intensity in terms of 
the Boltzmann value of the transport coefficient. 

4) The examples at the end of section 2 concern models satisfying the condi- 
tions of SDB, where the equilibrium distributions are completely factorized and 
position and velocity correlations are completely absent, as in (|ll|). However, 
this is not the case in LGCA's violating SDB, where computer simulations have 
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shown [Henon 1992, Bussemaker and Ernst 1992] the existence of (strong) on- 
node velocity correlations and (weaker) off-node spatial correlations in thermal 
equilibrium. The existence of such correlations is not yet understood. However, 
in driven diffusive systems [Zia et al 1993, Garrido et al 1990], where detailed 
balance is broken by imposing an external bias field or by boundary conditions, 
the absence of detailed balance does indeed give rise to long spatial correlations. 

5) As long as the analytic structure of the equilibrium correlation function 
g(x,y) is unknown, the fluctuation-dissipation relation ( fl0| ) is still valid but 
does not provide any information about the noise strength, B(x,y). The same 
remarks apply when dealing with stationary non-equilibrium states, support- 
ing stationary temperature gradients or shear rates. For real fluids the equal 
time correlation functions in such states are long ranged, and have been cal- 
culated from both Langevin and kinetic theories [Schmitz 1988] for which the 
noise strength is specified. As this information is not yet available for LGCA's 
further exploitation of ([ll]) is not possible. Similar limitations apply to (R) that 
describes fluctuations around non-stationary states. To calculate the correla- 
tion functions it is necessary to know the noise strength as a function of the 
non-equlibrium state, f(t). Again, this information is available for real fluids 
from kinetic theory but the corresponding analysis for LGCA's is incomplete 
[Ernst and Dufty 1994]. 

6) So far we have discussed stable systems. For unstable systems (phase sep- 
aration, spinodal deocmposition, pattern formation) the fundamental quantity 
is the equal time correlation function g(x, y, t) or its Fourier transform, the time 
dependent structure factor S(k,t). In the Cahn-Hilliard-Cook theory [Langer 
1992] the Langevin equation has been used to study the onset of instability, the 
wavelength of maximum growth, and initial patterns. Such studies have been 
done for reals systems as well as for LGCA's, quenched into a spatially uniform 
but metastable or unstable state. For some LGCA's simulations and analytic 
results (see Alexander et al [1992], Bussemaker and Ernst [1993], Ernst and 
Bussemaker 1994]) on S(k,t) and g(x,y,t) are becoming available, but studies 
on fluctuations around such states have yet to be done. Once B(x,y | f(t)) is 
understood on the basis of a microscopic model, that information may be used 
for mathematical modeling of g(x, y, t) via the fluctuation-dissipation relation 
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